clear; close all; clc;

%% Load GDP data

GDPdata = xlsread('Data_MATLAB.xlsx','GDP_PC','D1:BS15');

% Time series of the log of GDP per capita for each country
year =     GDPdata(1,:)'; 
ARG  = log(GDPdata(2,:))';
BRA  = log(GDPdata(3,:))';
ITA  = log(GDPdata(4,:))';
POR  = log(GDPdata(5,:))';
SPA  = log(GDPdata(6,:))';
% One-year differences
dyARG  = ARG(2:end) - ARG(1:end-1);
dyBRA  = BRA(2:end) - BRA(1:end-1);
dyITA  = ITA(2:end) - ITA(1:end-1);
dyPOR  = POR(2:end) - POR(1:end-1);
dySPA  = SPA(2:end) - SPA(1:end-1);
dy     = [dyITA,dyPOR,dySPA,dyARG,dyBRA];
year   = year(2:end);
ncountries = 5;

%% Set-ups for estimation

% Output: tables with mean, 5th and 95th percentiles of each parameter for
% each country
Ntheta     = 6;                             % number of parameters
table_mean = zeros(ncountries,Ntheta);      % means
table_p05  = zeros(ncountries,Ntheta);      % 5th percentiles
table_p95  = zeros(ncountries,Ntheta);      % 95th percentiles
NS = 2;  % number of states for Markov process
Nx = 2;  % length of unobserved vector
NC = 1;  % number of countries in each estimation (we estimate each country separately)
cnames = {'Italy','Portugal','Spain','Argentina','Brazil'}; 
% The scale parameters target rejection rates between 60% and 70%
scale_parameter = zeros(ncountries,1);
scale_parameter(1) = 0.005;      % Italy
scale_parameter(2) = 0.030;      % Portugal
scale_parameter(3) = 0.01358;    % Spain
scale_parameter(4) = 0.005;      % Argentina
scale_parameter(5) = 0.0065;     % Brazil
MH.CH    = 1;                    % number of chains to simulate
MH.Nburn = 25000;                % number of iterations for burn-in sample
MH.Nsim  = 225000;               % number of iterations to simulate draws
Nburn2 = 5000;                   % number of iterations to burn when computing the moments of the posterior distribution
MH.NB    = MH.Nburn + MH.Nsim;   % number of total iterations
MH.ext   = 10;                   % save a draw every "ext" elements of the chain
MH.print = 1;                    % print results if = 1 
MH.show  = 1000;
MH.save  = 1000;                

% Defining prior distributions - uniform
prior.p11lb = 0.05;    prior.p11ub = 0.99;  % bounds on the prior distribution for p11 (persistence of gL)
prior.p22lb = 0.05;    prior.p22ub = 0.99;  % bounds on the prior distribution for p22 (persistence of gH)
prior.mu1lb = -0.1;    prior.mu1ub = 0.1;   % bounds on the prior distribution for mu1 (gL)
prior.Dmulb =  0.0;    prior.Dmuub = 0.1;   % bounds on the prior distribution for m2-mu1; (gH = gL + (mu2-mu1))
prior.sg1lb =  0.001;  prior.sg1ub  = 0.075;% bounds on the prior distribution for sg (standard deviation of transitory shock)
prior.sg2lb =  0.001;  prior.sg2ub  = 0.075;% bounds on the prior distribution for sg (standard deviation of transitory shock)

%% Loop over countries

for iic = 1:ncountries
  
rng(1)
ydata = dy(:,iic);
cname  = cnames(iic);

% We select the initial point from the solution of the maximum likelihood estimation (MLE)
mu100 = mean(ydata(ydata<mean(ydata)));  % initial guess for gL to compute the MLE
mu200 = mean(ydata(ydata>=mean(ydata))); % initial guess for gH to compute the MLE
Dmu00 = mu200 - mu100;                   % initial guess for (gH-gL)
esim = ydata - (ydata<mean(ydata))*mu100 - (ydata>=mean(ydata))*mu200;
sg00 = std(esim)/2;                      % initial guess for sg (standard deviation of transitory shock) to compute the MLE
N11   = sum( (ydata(1:end-1)<mean(ydata)).*(ydata(2:end)<mean(ydata)) );
N1    = sum( (ydata<mean(ydata)) );
p1100 = N11/N1;                          % initial guess for p11 (persistence of gL) to compute the MLE
N22   = sum( (ydata(1:end-1)>=mean(ydata)).*(ydata(2:end)>=mean(ydata)) );
N2    = sum( (ydata>=mean(ydata)) );
p2200 = N22/N2;                          % initial guess for p22 (persistence of gH) to compute the MLE
% in the maximization/minimization we choose the weight between the lower bounds and upper bounds in the prior distribution of each parameter
% the weights must be between 0 and 1, so we choose log(weight/(1-weight));
% vector of weights in our initial guess:
htheta = [ (p1100-prior.p11ub)/(prior.p11lb-prior.p11ub); ...
           (p2200-prior.p22ub)/(prior.p22lb-prior.p22ub); ...
           (mu100-prior.mu1ub)/(prior.mu1lb-prior.mu1ub); ...
           (Dmu00-prior.Dmuub)/(prior.Dmulb-prior.Dmuub); ...
           (sg00 -prior.sg1ub )/(prior.sg1lb -prior.sg2ub ); ...
           (sg00 -prior.sg2ub )/(prior.sg1lb -prior.sg2ub )]; 
% initial guess for log(weight/(1-weight))
xx00    = log(htheta./(1-htheta));      
% Maximum Likelihood Estimation (MLE)
tr = 1; setmin = 1;
%options = optimset('Display','iter');
[xx_out,fval] = fminsearch(@(xx)posterior_fun_MA(xx,ydata,prior,NS,Nx,NC,tr,setmin),xx00);
% xx_out is our initial guess for log(weights/(1-weights)) for each
% parameter

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computing the MH Algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xx0      = xx_out;              % initial guess from the MLE
scale = scale_parameter(iic);   % scale parameter targeting rejection rate between 60% and 70%. 
% Covariance Matrix for draws --> update after burn-in periods:
SIGMA_xx    = scale*diag(max(abs(xx0),0.1));    
SIGMA_xx    = 0.5*(SIGMA_xx + SIGMA_xx');
xx_chain    = 189*ones(Ntheta,MH.NB);
lnL_chain   = 189*ones(MH.NB,1);
RAT_store   = 189*ones(MH.NB,1);
uu_draw     = unifrnd(0,1,MH.NB,1);  % for the acceptance/rejection decision
tr          = 1;   % transformation for posterior fun
setmin      = 0;   % = 1 if -PS in posterior fun  

% Initialize Chain
ib = 1;
xx_chain(:,ib) = xx0;
lnLch          = posterior_fun_MA(xx0,ydata,prior,NS,Nx,NC,tr,setmin);
lnL_chain(ib)  = lnLch;

% Markov Chain simulation
acc = 0; rej = 0; ibtot = 0; ss = MH.save; show = MH.show;
ch = 1; 
for ib = 2:MH.NB
    
    ibtot = ibtot + 1;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%    New draw Chain   %%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%   
    lnLj      = lnL_chain(ib-1);
    xx_new    = mvnrnd(xx_chain(:,ib-1),SIGMA_xx);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%   Update the Chain   %%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    lnL_new       = posterior_fun_MA(xx_new,ydata,prior,NS,Nx,NC,tr,setmin);                
    RAT_new       = exp(lnL_new-lnLj);         % exp(lnL_new)/exp(lnLj);        
    RAT_store(ib) = RAT_new;    
    if RAT_new > uu_draw(ib)  % this is a good theta draw -->  update
        xx_chain(:,ib) = xx_new;
        lnL_chain(ib)  = lnL_new;        
        acc = acc + 1;
    else
        xx_chain(:,ib) = xx_chain(:,ib-1);
        lnL_chain(ib)  = lnL_chain(ib-1);        
        rej = rej + 1;
    end
    
    if MH.print==1 && ch == 1   
        if show == 0
        disp(['MH   ',cell2mat(cnames(iic)),'  ib = ', num2str(ib),', acc = ', num2str(acc/(acc+rej),3),', rej = ', num2str(rej/(acc+rej),3), ', lnL = ', num2str(lnL_chain(ib),3),', eff = ', num2str((acc+rej)/ibtot,3)])
        show = MH.show;
        end
    end
    
    if ib == MH.Nburn  % reset acc and rej count
        disp('************************************************')
        disp('*****  resetting acc, rej and ibtot count  *****')
        disp('************************************************')
        acc = 0; rej = 0; ibtot = 0;
    end
    
    show = show - 1;
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Posterior computations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% transform xx into theta
theta_chain = xx_chain;
alphap11 = exp(xx_chain(1,:))./(1+exp(xx_chain(1,:))); theta_chain(1,:)   = alphap11.*prior.p11lb + (1-alphap11).*prior.p11ub; % p1
alphap22 = exp(xx_chain(2,:))./(1+exp(xx_chain(2,:))); theta_chain(2,:)   = alphap22.*prior.p22lb + (1-alphap22).*prior.p22ub; % p2
alphamu1 = exp(xx_chain(3,:))./(1+exp(xx_chain(3,:))); theta_chain(3,:)   = alphamu1.*prior.mu1lb + (1-alphamu1).*prior.mu1ub; % gH
alphaDmu = exp(xx_chain(4,:))./(1+exp(xx_chain(4,:))); theta_chain(4,:)   = alphaDmu.*prior.Dmulb + (1-alphaDmu).*prior.Dmuub; % gH-gL
alphasg1 = exp(xx_chain(5,:))./(1+exp(xx_chain(5,:))); theta_chain(5,:)   = alphasg1.*prior.sg1lb + (1-alphasg1).*prior.sg1ub; % sg1
alphasg2 = exp(xx_chain(6,:))./(1+exp(xx_chain(6,:))); theta_chain(6,:)   = alphasg2.*prior.sg2lb + (1-alphasg2).*prior.sg2ub; % sg2
Dmu_chain        = theta_chain(4,:);           
theta_chain(4,:) = theta_chain(3,:) + Dmu_chain;  % transforming (gH-gL) into gH
% save entire chain
theta_chain_lng = theta_chain;
% discard MH.Nburn initial draws, and then keep one every MH.ext draws
in          = MH.Nburn+1:MH.ext:(MH.Nburn+MH.Nsim);
theta_chain = theta_chain_lng(:,in);

MH.Nfn = size(theta_chain,2);  
count  = 1:1:MH.Nfn;
recmeans = 189*ones(Ntheta,MH.Nfn);
for ii = 1:Ntheta
    recmeans(ii,:) = cumsum(theta_chain(ii,:))./count;
end


% Computing mean and 5th and 95th percentiles
table_mean(iic,:) = mean(theta_chain(:,Nburn2:end),2)';
table_p05(iic,:)  = prctile(theta_chain(:,Nburn2:end)',5);
table_p95(iic,:)  = prctile(theta_chain(:,Nburn2:end)',95);

figure
subplot(2,3,1)
hist(theta_chain(1,:),100)
title(['$\hat{p}_{11} = $', num2str(mean(theta_chain(1,:)),3)],'Interpreter','LaTex','Fontsize',32)
xlabel('$p_{11}$','Interpreter','LaTex')
set(gca,'Xgrid','on','Ygrid','on','Fontsize',21)
subplot(2,3,4)
hist(theta_chain(2,:),100)
title(['$\hat{p}_{22} = $', num2str(mean(theta_chain(2,:)),3)],'Interpreter','LaTex','Fontsize',32)
xlabel('$p_{22}$','Interpreter','LaTex')
set(gca,'Xgrid','on','Ygrid','on','Fontsize',21)
subplot(2,3,2)
hist(theta_chain(3,:),100)
title(['$\hat{\mu}_1 = $', num2str(mean(theta_chain(3,:)),3)],'Interpreter','LaTex','Fontsize',32)
xlabel('$\mu_1$','Interpreter','LaTex')
set(gca,'Xgrid','on','Ygrid','on','Fontsize',21)
subplot(2,3,5)
hist(theta_chain(4,:),100)
title(['$\hat{\mu}_2 = $', num2str(mean(theta_chain(4,:)),3)],'Interpreter','LaTex','Fontsize',32)
xlabel('$\mu_2$','Interpreter','LaTex')
set(gca,'Xgrid','on','Ygrid','on','Fontsize',21)
subplot(2,3,3)
hist(theta_chain(5,:),100)
title(['$\hat{\sigma_1} = $', num2str(mean(theta_chain(5,:)),3)],'Interpreter','LaTex','Fontsize',32)
xlabel('$\sigma$','Interpreter','LaTex')
set(gca,'Xgrid','on','Ygrid','on','Fontsize',21)
subplot(2,3,6)
hist(theta_chain(6,:),100)
title(['$\hat{\sigma_2} = $', num2str(mean(theta_chain(6,:)),3)],'Interpreter','LaTex','Fontsize',32)
xlabel('$\sigma$','Interpreter','LaTex')
set(gca,'Xgrid','on','Ygrid','on','Fontsize',21)

figure
subplot(2,3,1)
plot(recmeans(1,:),'r','LineWidth',3)
title(['Rec Mean $p_{11}$ -- $\hat{p}_{11} = $', num2str(mean(theta_chain(1,:)),3)],'Interpreter','LaTex','Fontsize',32)
xlabel('Iteration','Interpreter','LaTex')
xlim([0 MH.Nfn])
set(gca,'Xgrid','on','Ygrid','on','Fontsize',21)
subplot(2,3,4)
plot(recmeans(2,:),'r','LineWidth',3)
title(['Rec Mean $p_{22}$ -- $\hat{p}_{22} = $', num2str(mean(theta_chain(2,:)),3)],'Interpreter','LaTex','Fontsize',32)
xlabel('Iteration','Interpreter','LaTex')
xlim([0 MH.Nfn])
set(gca,'Xgrid','on','Ygrid','on','Fontsize',21)
subplot(2,3,2)
plot(recmeans(3,:),'r','LineWidth',3)
title(['Rec Mean $\mu_1$ -- $\hat{\mu}_1 = $', num2str(mean(theta_chain(3,:)),3)],'Interpreter','LaTex','Fontsize',32)
xlabel('Iteration','Interpreter','LaTex')
xlim([0 MH.Nfn])
set(gca,'Xgrid','on','Ygrid','on','Fontsize',21)
subplot(2,3,5)
plot(recmeans(4,:),'r','LineWidth',3)
title(['Rec Mean $\mu_2$ -- $\hat{\mu}_2 = $', num2str(mean(theta_chain(4,:)),3)],'Interpreter','LaTex','Fontsize',32)
xlabel('Iteration','Interpreter','LaTex')
xlim([0 MH.Nfn])
set(gca,'Xgrid','on','Ygrid','on','Fontsize',21)
subplot(2,3,3)
plot(recmeans(5,:),'r','LineWidth',3)
title(['Rec Mean $\sigma_1$ -- $\hat{\sigma_1} = $', num2str(mean(theta_chain(5,:)),3)],'Interpreter','LaTex','Fontsize',32)
xlabel('Iteration','Interpreter','LaTex')
xlim([0 MH.Nfn])
set(gca,'Xgrid','on','Ygrid','on','Fontsize',21)
subplot(2,3,6)
plot(recmeans(6,:),'r','LineWidth',3)
title(['Rec Mean $\sigma_2$ -- $\hat{\sigma_2} = $', num2str(mean(theta_chain(6,:)),3)],'Interpreter','LaTex','Fontsize',32)
xlabel('Iteration','Interpreter','LaTex')
xlim([0 MH.Nfn])
set(gca,'Xgrid','on','Ygrid','on','Fontsize',21)

end

filename = 'Table_A1.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)

% Table A1
fprintf('----------------------------------------------------------------------- \n')
fprintf('TABLE A1: Posterior distributions: different σ across growth states \n')
fprintf('----------------------------------------------------------------------- \n')
fprintf('           \t   gL   \t  gH   \t   pL \t  pH   \t  sgL \t  sgH   \n')
for iic=1:ncountries
fprintf([cell2mat(cnames(iic)),'\n'])
fprintf(' mean     \t  %5.3f \t %5.3f \t %5.3f \t %5.3f \t %5.3f \t %5.3f \n', [table_mean(iic,3),table_mean(iic,4),table_mean(iic,1),table_mean(iic,2),table_mean(iic,5),table_mean(iic,6)])
fprintf(' p05     \t  %5.3f \t %5.3f \t %5.3f \t %5.3f \t %5.3f \t %5.3f \n', [table_p05(iic,3),table_p05(iic,4),table_p05(iic,1),table_p05(iic,2),table_p05(iic,5),table_p05(iic,6)])
fprintf(' p95     \t  %5.3f \t %5.3f \t %5.3f \t %5.3f \t %5.3f \t %5.3f \n', [table_p95(iic,3),table_p95(iic,4),table_p95(iic,1),table_p95(iic,2),table_p95(iic,5),table_p95(iic,6)])
fprintf('\n')
end
fprintf('----------------------------------------------------------------------- \n')

diary off
        
